library(knitr)
opts_chunk$set(warning = FALSE, message = FALSE, fig.path = 'figs/', dev.args = list(bg = 'transparent'), eval = T)

library(raster)
library(sf)
library(tidyverse)
library(mapview)
library(proj4shortcut)
library(plotly)

prj_geo <- geo_wgs84
prj_pro <- utm_wgs84('11s')

August 21st

Map of in situ samples with chlorophyll concentration.

# get coverage of in situ points by date
insit <- read_csv('raw/CHL_2017_08_21.csv', col_names = F) %>% 
  rename(
    site = X1, 
    lat = X2, 
    lon = X3, 
    chl = X4
  ) %>% 
  mutate(lat = ifelse(site == 34, 33.66954, lat)) 
coordinates(insit) <- ~ lon + lat
proj4string(insit) <- prj_geo
mapview(insit, zcol = 'chl', legend = T) 

Overlap of in situ samples (red numbers) with the spatial extent of each image. Note that the spatial extent is rectangular and not all pixels in the extent have values.

# transform insit to projected
insit <- spTransform(insit, CRS(prj_pro)) 

# raster tiffs
imgs_pth <- list.files('img/GeoTiff_08_21/', pattern = '^Raster*', full.names = T)

# get raster extents
out <- list()
cmbs <- for(i in 1:length(imgs_pth)){
  
  # raster aggregate, dissolve to polygon
  rst_chk <- stack(imgs_pth[[i]]) %>% 
    .[[1]] %>% 
    aggregate(fact = 10, fun = mean)
  values(rst_chk) <- ifelse(is.na(values(rst_chk)), NA, 1)
  ply <- rst_chk %>% 
    rasterToPolygons(dissolve = T)
  
  out[[i]] <- ply
  
}

# fortify polygons for plot
allply <- out %>% 
  map(., function(x){
    dat <- fortify(x) %>% 
      mutate(fl = gsub('\\.[0-9]$', '', names(x@data)))
    return(dat)
    }) %>% 
  enframe %>% 
  unnest

# plot
p <- ggplot() + 
  geom_polygon(data = allply, aes(x = long, y = lat, group = fl), alpha = 0.5) + 
  geom_text(data = data.frame(insit), aes(x = lon, y = lat, label = site)) + 
  coord_equal()
ggplotly(p)

Extracted pixel values for each of four bands in each image where sample sites occurred within an image.

# extract values from raster
out <- list()
for(i in 1:length(imgs_pth)){
  
  # cat(i, 'of', length(imgs_pth), '\n')
  
  # extract raster cells by points in insit
  rst <- stack(imgs_pth[[i]])
  crs(rst) <- prj_pro
  rst_chk  <- rst %>% 
    raster::extract(insit, buffer = 0) %>% 
    enframe %>% 
    bind_cols(insit@data, .) %>% 
    select(-name) %>% 
    filter(map(.$value, function(x) !anyNA(x)) %>% unlist) %>% 
    mutate(value = map(value, function(x){

      sumpt <- x %>% 
        data.frame(din = .) %>% 
        rownames_to_column('img_bnd') %>% 
        filter(!grepl('\\.5$', img_bnd))
      return(sumpt)
      
    })) %>% 
    unnest

  # append to output
  out[[i]] <- rst_chk

}

# combine all extracted samples
# get ndvi, gndvi, vigreen
# current not in nm
exts <- out %>% 
  do.call('rbind', .) %>% 
  remove_rownames() %>% 
  separate(img_bnd, c('img', 'bnd'), sep = '\\.') %>% 
  select(-img) %>% 
  mutate(bnd = factor(bnd, levels = c('1', '2', '3', '4'), labels = c('grn', 'red', 'edg', 'nir'))) %>% 
  spread(bnd, din) %>% 
  mutate(
    ndvi = (nir - red) / (nir + red), # normalized diff veg index
    gndvi = (nir - grn) / (nir + grn), # green normalized diff veg index
    vigrn = (grn - red) / (grn + red) # normalized diff of green and red bands
  ) %>% 
  gather('var', 'val', -site, -chl)
exts
## # A tibble: 91 x 4
##     site   chl var     val
##    <int> <dbl> <chr> <dbl>
##  1     1 146   grn   0.267
##  2     2 143   grn   0.271
##  3    16  52.6 grn   0.269
##  4    17 209   grn   0.308
##  5    31 189   grn   0.330
##  6    32 157   grn   0.315
##  7    33 185   grn   0.337
##  8    34 170   grn   0.304
##  9    35 186   grn   0.302
## 10    36 285   grn   0.647
## # ... with 81 more rows

Scatterplots of band values (din) with measured chlorophyll.

ggplot(exts, aes(x = val, y = chl)) + 
  geom_point() + 
  facet_wrap(~var, ncol = 2, scales = 'free_x') + 
  stat_smooth(method = 'lm') + 
  theme_bw()

Working with model selection:

library(caret)
library(MuMIn)
# https://pix4d.com/sequoia-faq/#3

tomod <- exts %>% 
  mutate(bnd = paste0('b', bnd)) %>% 
  spread(bnd, din)

glb <- lm(chl ~ (b1 + b2 + b2 + b3 + b4)^3 + I(b1^2) + I(b2^2) + I(b3^2) + I(b4^2), tomod, 
          na.action = 'na.pass')

tmp <- dredge(glb, evaluate = F)

# createFolds(1:nrow(tomod), k = 5)

September 6th

Map of in situ samples with chlorophyll concentration.

# get coverage of in situ points by date
insit <- read_csv('raw/CHL_2017_09_06.csv', col_names = F) %>% 
  rename(
    site = X1, 
    lat = X2, 
    lon = X3, 
    chl = X4
  ) %>% 
  mutate(lon = ifelse(site == 39, -117.3597, lon))
coordinates(insit) <- ~ lon + lat
proj4string(insit) <- prj_geo
mapview(insit, zcol = 'chl', legend = T) 

Overlap of in situ samples (red numbers) with the spatial extent of each image. Note that the spatial extent is rectangular and not all pixels in the extent have values.

# transform insit to projected
insit <- spTransform(insit, CRS(prj_pro)) 

# raster tiffs
imgs_pth <- list.files('img/GeoTiff_09_06/', pattern = '*\\.tif$', full.names = T)

# get raster extents
out <- list()
cmbs <- for(i in 1:length(imgs_pth)){
  
  # raster aggregate, dissolve to polygon
  rst_chk <- stack(imgs_pth[[i]]) %>% 
    .[[1]] %>% 
    aggregate(fact = 10, fun = mean)
  values(rst_chk) <- ifelse(is.na(values(rst_chk)), NA, 1)
  ply <- rst_chk %>% 
    rasterToPolygons(dissolve = T)
  
  out[[i]] <- ply
  
}

# fortify polygons for plot
allply <- out %>% 
  map(., function(x){
    dat <- fortify(x) %>% 
      mutate(fl = gsub('\\.[0-9]$', '', names(x@data)))
    return(dat)
    }) %>% 
  enframe %>% 
  unnest

# plot
p <- ggplot() + 
  geom_polygon(data = allply, aes(x = long, y = lat, group = fl), alpha = 0.5) + 
  geom_text(data = data.frame(insit), aes(x = lon, y = lat, label = site)) + 
  coord_equal()
ggplotly(p)

Extracted pixel values for each of four bands in each image where sample sites occurred within an image.

# extract values from raster
out <- list()
for(i in 1:length(imgs_pth)){
  
  # cat(i, 'of', length(imgs_pth), '\n')
  
  rst <- stack(imgs_pth[[i]])
  crs(rst) <- prj_pro
  rst_chk  <- rst %>% 
    raster::extract(insit) %>% 
    data.frame(insit@data, .) %>%
    gather('img_bnd', 'din', -site, -chl) %>% 
    na.omit

  # out <- c(out, list(tmp))
  out[[i]] <- rst_chk

}

exts <- out %>% 
  do.call('rbind', .) %>% 
  remove_rownames() %>% 
  separate(img_bnd, c('img', 'bnd'), sep = '\\.') %>% 
  group_by(site, chl, bnd) %>% 
  summarise(din = mean(din)) %>% 
  ungroup %>% 
  mutate(bnd = factor(bnd, levels = c('1', '2', '3', '4'), labels = c('grn', 'red', 'edg', 'nir'))) %>% 
  spread(bnd, din) %>% 
  mutate(
    ndvi = (nir - red) / (nir + red), # normalized diff veg index
    gndvi = (nir - grn) / (nir + grn), # green normalized diff veg index
    vigrn = (grn - red) / (grn + red) # normalized diff of green and red bands
  ) %>% 
  gather('var', 'val', -site, -chl)
exts
## # A tibble: 70 x 4
##     site   chl var     val
##    <int> <dbl> <chr> <dbl>
##  1     1  62.6 grn   0.319
##  2     2 487   grn   0.322
##  3     3 185   grn   0.315
##  4    16 380   grn   0.326
##  5    17 295   grn   0.298
##  6    18 306   grn   0.305
##  7    19 293   grn   0.313
##  8    31 254   grn   0.384
##  9    32 215   grn   0.338
## 10    33 152   grn   0.497
## # ... with 60 more rows

Scatterplots of band values (din) with measured chlorophyll.

ggplot(exts, aes(x = val, y = chl)) + 
  geom_point() + 
  facet_wrap(~var, ncol = 2, scales = 'free_x') + 
  stat_smooth(method = 'lm') + 
  theme_bw()